Jupyter Notebook for Julia¶

A. Install Julia in your Google Colab.¶

If Julia is already isntalled in your Google Colab, go to B. If Julia and Jupyter Notebook are installed in your computer, go to D.

"Run" means

  • Put the cursor on the cel and shift+enter on Jupyter Notebook
  • Tap the circle near the upper left corner of the cell on Google Colab

Run the following cell.¶

In [ ]:
!curl -fsSL https://install.julialang.org | sh -s -- --yes
!julia -v

B. Choose Julia in Runtime.¶

Go to "Runtime" $\rightarrow$ "Change runtime type", and choose "Julia".¶

C. Install packages in your Julia.¶

If the following packages are already installed, go to D.

Run the following cell. The installation takes 5-10 minutes.¶

In [ ]:
using Pkg
Pkg.add("LaTeXStrings")
Pkg.add("Images")
Pkg.add("Plots")
Pkg.add("Colors")
Pkg.add("Downloads")
Pkg.add("FileIO")
Pkg.add("LowRankApprox")
Pkg.add("Wavelets")
Pkg.add("DSP")
Pkg.add("Dates")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("Statistics")

D. Load packages.¶

Run the following cell.¶

In [2]:
using Random
using LinearAlgebra
using LaTeXStrings
using Images
using Plots
using Colors
using Downloads
using FileIO
using LowRankApprox
using Wavelets
using DSP
using Dates
using CSV
using DataFrames
using Statistics
In [ ]:
# space

Topics¶

  1. The Collatz conjecture
  2. RGB images
  3. RGB movies
  4. Singular value decomposition and low rank approximation
  5. Wavelet decomposition
  6. Public data and visualization
  7. Central limit theorem for the binomial distribution $B(n,p)$
  8. Taylor series of $\sin{x}$
  9. Fourier series of a step function and a sawtooth function
  10. Weyl's equidistribution theorem

1. The Collatz conjecture¶

Generate a sequence $\{x[n]\}_{n=1}^\infty$ by the following process:

  1. Given an arbitrary positive integer $m$, set $x[1]:=m$.
  2. For $n=2,3,4,\dotsc$, if $x[n-1]:=1$, stop the process, otherwise

\begin{aligned} x[n] & :=\dfrac{x[n-1]}{2}\quad (x[n-1]\ \text{is even}), \\ x[n] & := 3x[n-1]+1\quad (x[n-1]\ \text{is odd}). \end{aligned}

In 1937 Lothar Collatz introduce the following conjecture:

For any $m=1,2,3,\dotsc$, there exists $n_m=1,2,3,\dotsc$ such that $x[n_m]=1$.¶

See this paper for instance. Here are some examples:

$1, \quad 2\rightarrow1, \quad 3\rightarrow10\rightarrow5\rightarrow16\rightarrow8\rightarrow4\rightarrow2\rightarrow1, \quad 4\rightarrow2\rightarrow1, \quad\dotsc$

Barina verified the Collatz conjecture is true for all $m\leqq2^{68}$ using workstations. See this paper.

We examine this conjecture for $m=1,2,3,\dotsc,500$.

Run the following cell to create Collatz sequences for $x[1]=m=1,\dotsc,500$, and observe the behavior of each sequence.¶

In [27]:
# upper bound of m
M=500;

x=ones(Int64,150);

animcollatz = @animate for m=1:M

x[1]=m
for n=2:150
    if x[n-1]==1
       x[n]=x[n]
    else
        if mod(x[n-1],2)==0
            x[n]=Int64.(x[n-1]/2)
        else x[n]=3*x[n-1]+1
        end
    end
end
	
scatter(x,
	 ylim=(1,2000),
	 title="Collatz sequence x[n] with x[1]=$m",
	 xlabel=L"n",
	 ylabel=L"x[n]",
	 xticks=([1,20,40,60,80,100,120,140]),
	 legend=false,
     clear=true)

end

gif(animcollatz, "collatz_animation.gif", fps = 5)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/collatz_animation.gif
Out[27]:
No description has been provided for this image

Run the following cell to plot the pairs of $(m,n_m)$ for $m=1,\dotsc,500$ in a scatter graph.¶

In [24]:
function collatz(n)
    n % 2 == 0 ? n ÷ 2 : 3n + 1
end

function stopping_time(m)
        x, steps = m, 0
        while x > 1
              x = collatz(x)
              steps += 1
		end
			return steps
end

N = 1:500
steps = [stopping_time(m) for m in N]

scatter(N, steps, 
		legend=false,
		xlabel=L"\mathrm{Initial}\ \mathrm{value}\ \ m", 
		ylabel=L"\mathrm{Steps}\ \ n_m\ \  \mathrm{to}\ \ 1",
        title="Collatz stopping steps", markersize=2)
Out[24]:
No description has been provided for this image

2. RGB images¶

The data of an RGB image is a triplet of matrices of the same size. This is essentially a triplet of grayscale images of the same size. These three grayscale images, colored red, green, and blue, are then overlaid to form an RGB image.

Run the following cell to create three $m{\times}n=10\times15$ matrices $R,G,B$ by random numbers.¶

In [5]:
Random.seed!(1)
m, n = 10, 15 # size of matrices

Z = zeros(m, n); # zero matrix
R = rand(m, n);
G = rand(m, n);
B = rand(m, n) # See below
Out[5]:
10×15 Matrix{Float64}:
 0.670414   0.643352  0.376138   0.983634   …  0.899505  0.380349   0.221796
 0.691898   0.100056  0.218766   0.40892       0.424926  0.958143   0.950771
 0.192718   0.164962  0.60248    0.0424498     0.237027  0.610986   0.839523
 0.516169   0.959669  0.753166   0.705895      0.189413  0.619391   0.736151
 0.0475605  0.72755   0.852577   0.0647826     0.049791  0.849654   0.777125
 0.245685   0.436295  0.88632    0.792918   …  0.587347  0.515305   0.0444534
 0.979139   0.396001  0.0642346  0.080433      0.171102  0.0882561  0.0745833
 0.0294498  0.419402  0.890999   0.998768      0.265295  0.811548   0.626937
 0.224011   0.525583  0.988467   0.0620525     0.480608  0.0926662  0.601864
 0.80577    0.21062   0.292272   0.199147      0.582327  0.494502   0.440245

Run the following cell to colorize $R,G,B$, and show red, green, blue and RGB images.¶

In [6]:
imgR = colorview(RGB, R, Z, Z)
imgG = colorview(RGB, Z, G, Z)
imgB = colorview(RGB, Z, Z, B)
imgRGB = colorview(RGB, R, G, B)

p1 = plot(imgR, grid=false, axis=false, title="Red")
p2 = plot(imgG, grid=false, axis=false, title="Green")
p3 = plot(imgB, grid=false, axis=false, title="Blue")
p4 = plot(imgRGB, grid=false, axis=false, title="RGB")

plot(p1, p2, p3, p4, layout=(2,2), size=(600,450))
Out[6]:
No description has been provided for this image

3. RGB movies¶

The data of an RGB movie is the finite number of triplets of matrices of the same size. This is essentially a sequence of RGB images of the same size.

Run the following cell to create 100 $m{\times}n=10\times15$ RGB images by random numbers.¶

In [7]:
Random.seed!(1)
m, n = 10, 15 # size of matrices
K = 100 # steps, length of the sequence

Z = zeros(m, n);
RR = rand(m, n, K);
GG = rand(m, n, K);
BB = rand(m, n, K);

Run the following cell to see the movie.¶

In [8]:
animrgb = @animate for k=1:K

movieR = colorview(RGB, RR[:,:,k], Z, Z)
movieG = colorview(RGB, Z, GG[:,:,k], Z)
movieB = colorview(RGB, Z, Z, BB[:,:,k])
movieRGB = colorview(RGB, RR[:,:,k], GG[:,:,k], BB[:,:,k])

p1 = plot(movieR, grid=false, axis=false, title="Red")
p2 = plot(movieG, grid=false, axis=false, title="Green")
p3 = plot(movieB, grid=false, axis=false, title="Blue")
p4 = plot(movieRGB, grid=false, axis=false, title="RGB")

plot(p1, p2, p3, p4, layout=(2,2), size=(600,450))

end

gif(animrgb, "rgb_animation.gif", fps = 3)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/rgb_animation.gif
Out[8]:
No description has been provided for this image

4. Singular value decomposition and low rank approximation¶

In general an $m{\times}n$ matrix $A=\begin{bmatrix}\vec{a}_1 & \dotsb & \vec{a}_n \end{bmatrix}$ has a specific nonnegative integer which is called as the rank of $A$ defined by

$$r := \operatorname{dim}\bigl( \{ A\vec{x} = x_1\vec{a}_1+\dotsb+x_n\vec{a}_n \ \vert \ \vec{x} \in \mathbb{R}^n \} \bigr).$$

It follows that $1 \leqq r \leqq \min\{m,n\}$ unless all the entries of $A$ are zero. It is known that there exist

$$\sigma_1 \geqq \dotsb \geqq \sigma_r>0, \quad \vec{u}_1,\dotsc,\vec{u}_r \in \mathbb{R}^m, \quad \vec{v}_1,\dotsc,\vec{v}_r \in \mathbb{R}^n$$

such that we have the singular value decomposition:

$$A=\sum_{j=1}^r\sigma_j\vec{u}_j\vec{v}_j^T.$$

Set

$$A_k:=\sum_{j=1}^k\sigma_j\vec{u}_j\vec{v}_j^T, \quad k=1,\dotsc,r.$$

It is known that the rank of $A_k$ is $k$ and closest to $A$ among all the $m{\times}n$ matrices of rank $k$, that is, $A_k$ is the best approximation of $A$ among all the $m{\times}n$ matrices of rank $k$.

We compute a $432\times768\times3$ RGB image of hibiscus. All the R, G, B matrices are of rank $432$.

Run the following cell to compute the SVD of the RGB-elements of the picture.¶

In [9]:
url = "https://raw.githubusercontent.com/fiomfd/hands-on/refs/heads/main/data/hibiscus.jpg"
Downloads.download(url, "hibiscus.jpg")
I=load("./hibiscus.jpg");
X=imresize(I, ratio=1/5);
(p,q)=size(X);
AAA=channelview(X);
RRR=Array{Float64}(AAA[1,:,:]);
GGG=Array{Float64}(AAA[2,:,:]);
BBB=Array{Float64}(AAA[3,:,:]);

RU, RS, RV=psvd(RRR);
GU, GS, GV=psvd(GGG);
BU, BS, BV=psvd(BBB);
	
rank=60;
DR=zeros(p,q,rank);
DG=zeros(p,q,rank);
DB=zeros(p,q,rank);
for r=1:rank
    DR[:,:,r]=sum(RS[n]*RU[1:p,n]*(RV[1:q,n])' for n=1:r);
    DG[:,:,r]=sum(GS[n]*GU[1:p,n]*(GV[1:q,n])' for n=1:r);
    DB[:,:,r]=sum(BS[n]*BU[1:p,n]*(BV[1:q,n])' for n=1:r);
end

Run the following cell to see the movie of the best low rank approximation for $r=1,\dotsc,60$.¶

In [10]:
animsvd = @animate for k=1:rank

W=zeros(3,p,q,rank);
for r=1:rank
	W[1,:,:,r]=DR[:,:,r];
    W[2,:,:,r]=DG[:,:,r];
    W[3,:,:,r]=DB[:,:,r];
end

SVD1=plot(colorview(RGB,W[:,:,:,k]),
        title="Rank-$k approximation",
        titlefontsize=12,
        xaxis=false, 
        xticks=false, 
        yaxis=false, 
        yticks=false, 
        grid=false);
SVD2=plot(colorview(RGB,X),
        title="Original RGB image (rank 432)",
        titlefontsize=12,
        xaxis=false, 
        xticks=false, 
        yaxis=false, 
        yticks=false, 
        grid=false);
plot(SVD1,SVD2, layout=(1,2), size=(650,220))

end

gif(animsvd, "svd_animation.gif", fps = 3)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/svd_animation.gif
Out[10]:
No description has been provided for this image

5. Wavelet decomposition¶

The discrete wavelet is a pair of orthonomal two vectors $\vec{u}$ and $\vec{v}$ satisfying some condition, and the most typical example of the discrete wavelets is the Haar wavelet. They are very simple vectors:

$$\vec{u} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}, \quad \vec{v} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \in\mathbb{R}^{N},$$

Roughly speaking the discrete wavelet transform by the Haar wavelet is to take averages of neighboring elements. The elements of a vector are transformed to the average of neighboring two elements, and the elements of a matrix are converted to the average of neighboring 2 by 2 elements. The following examples show the outcomes of repeating the discere wavelet tranform twice:

$$\vec{a}_0 = \begin{bmatrix} 1 \\ 3 \\ 5 \\ 7 \end{bmatrix} \mapsto \vec{a}_1 = \begin{bmatrix} 2 \\ 2 \\ 6 \\ 6 \end{bmatrix} \mapsto \vec{a}_2 = \begin{bmatrix} 4 \\ 4 \\ 4 \\ 4 \end{bmatrix}$$

$$A_0 = \begin{bmatrix} 0 & 2 & 4 & 6 \\ 8 & 10 & 12 & 14 \\ 16 & 18 & 20 & 22 \\ 24 & 26 & 28 & 30 \end{bmatrix} \mapsto A_1 = \begin{bmatrix} 5 & 5 & 9 & 9 \\ 5 & 5 & 9 & 9 \\ 21 & 21 & 25 & 25 \\ 21 & 21 & 25 & 25 \end{bmatrix} \mapsto A_2 = \begin{bmatrix} 15 & 15 & 15 & 15 \\ 15 & 15 & 15 & 15 \\ 15 & 15 & 15 & 15 \\ 15 & 15 & 15 & 15 \end{bmatrix}.$$

The output such as $\vec{a}_\ell$ and $A_\ell$ of $\ell$ times operations is said to be the approximation part of level $\ell$, and the remainder term such as $\vec{a}_0-\vec{a}_\ell$ and $A_0-A_\ell$ is called the detail part of level $\ell$. In terms of Fourier analysis, the approximation part is the low frequency part, and the detail part is the high frequency part and is similar to the edge detection.

We observe the wavelet decomposition of the $384\times512$ RGB image of fried noodle using the Haar wavelet. We can get the wavelet decomposition up to level $7$ since $384=2^7\times3$ and $512=2^9$.

Run the following cell to implement the wavelet decomposition.¶

In [11]:
url = "https://raw.githubusercontent.com/fiomfd/hands-on/refs/heads/main/data/CityU.jpg"
Downloads.download(url, "CityU.jpg")
I8=load("./CityU.jpg");
X8=imresize(I8, ratio=1/8);
(p8,q8)=size(X8);
A8=channelview(X8);
R8=Array{Float64}(A8[1,:,:]);
G8=Array{Float64}(A8[2,:,:]);
B8=Array{Float64}(A8[3,:,:]);

# Maximum level
L=7;

# Decomposition Filter 
XR=zeros(p8,q8,L);
XG=zeros(p8,q8,L);
XB=zeros(p8,q8,L);
for l=1:L
	XR[:,:,l]=dwt(R8, wavelet(WT.haar), l);
	XG[:,:,l]=dwt(G8, wavelet(WT.haar), l);
	XB[:,:,l]=dwt(B8, wavelet(WT.haar), l);
end

# Splitting 
XRapprox=zeros(p8,q8,L);
XGapprox=zeros(p8,q8,L);
XBapprox=zeros(p8,q8,L);	
for l=1:L
	XRapprox[1:Int(p8/2^l),1:Int(q8/2^l),l]=XR[1:Int(p8/2^l),1:Int(q8/2^l),l];
    XGapprox[1:Int(p8/2^l),1:Int(q8/2^l),l]=XG[1:Int(p8/2^l),1:Int(q8/2^l),l];
	XBapprox[1:Int(p8/2^l),1:Int(q8/2^l),l]=XB[1:Int(p8/2^l),1:Int(q8/2^l),l];
end
XRdetail=XR-XRapprox;
XGdetail=XG-XGapprox;
XBdetail=XB-XBapprox;

# Composition Filter
YRapprox=zeros(p8,q8,L);
YGapprox=zeros(p8,q8,L);
YBapprox=zeros(p8,q8,L);
YRdetail=zeros(p8,q8,L);
YGdetail=zeros(p8,q8,L);
YBdetail=zeros(p8,q8,L);

for l=1:L
	YRapprox[:,:,l]=idwt(XRapprox[:,:,l], wavelet(WT.haar), l);
	YGapprox[:,:,l]=idwt(XGapprox[:,:,l], wavelet(WT.haar), l);
	YBapprox[:,:,l]=idwt(XBapprox[:,:,l], wavelet(WT.haar), l);
	YRdetail[:,:,l]=idwt(XRdetail[:,:,l], wavelet(WT.haar), l);
	YGdetail[:,:,l]=idwt(XGdetail[:,:,l], wavelet(WT.haar), l);
	YBdetail[:,:,l]=idwt(XBdetail[:,:,l], wavelet(WT.haar), l);
end

# RGB
Wapprox=zeros(3,p8,q8,L+1);
Wdetail=zeros(3,p8,q8,L+1);
Wapprox[:,:,:,1]=A8;
for l=2:L+1
	Wapprox[1,:,:,l]=YRapprox[:,:,l-1];
	Wapprox[2,:,:,l]=YGapprox[:,:,l-1];
	Wapprox[3,:,:,l]=YBapprox[:,:,l-1];
	Wdetail[1,:,:,l]=YRdetail[:,:,l-1];
	Wdetail[2,:,:,l]=YGdetail[:,:,l-1];
	Wdetail[3,:,:,l]=YBdetail[:,:,l-1];
end

Run the previous cell to see a movie of the wavelet decomposition.¶

In [12]:
animhaar = @animate for l=0:L

P1=plot(colorview(RGB,A8[:,:,:]),
    title="Original RGB image",
    titlefontsize=12,
    xaxis=false, 
    xticks=false, 
    yaxis=false, 
    yticks=false, 
    grid=false);
P2=plot(colorview(RGB,Wapprox[:,:,:,l+1]),
    title="Level $l Approximation",
    titlefontsize=12,
    xaxis=false, 
    xticks=false, 
    yaxis=false, 
    yticks=false, 
    grid=false);
P3=plot(colorview(RGB,Wdetail[:,:,:,l+1]),
    title="Level $l Detail",
    titlefontsize=12,
    xaxis=false, 
    xticks=false, 
    yaxis=false, 
    yticks=false, 
    grid=false);
plot(P1,P2,P3,size=(900,240),layout=(1,3))

end

gif(animhaar, "svd_animation.gif", fps = 0.2)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/svd_animation.gif
Out[12]:
No description has been provided for this image

6. Public data and visualization¶

Let $\vec{x}=\bigl[x[1],\dotsc,x[N]\bigl]^T\in\mathbb{R}^N$ be an $N$-dimensional real column vector. This can describes data such as $x[i]$ means the highest temperature of the $i$-th day for $i=1,\dotsc,N$. One can exploit some information such as the average $m$ and the variance $\sigma^2$:

$$m:=\frac{1}{N}\sum_{i=1}^Nx[i], \quad \sigma^2:=\frac{1}{N}\sum_{i=1}^N(x[i]-m)^2.$$

In what follows Julia downloads open data, and visualizes information taking averages, maximums and minimums. More concretely Julia downloads

  • the daily maximum tenperature
  • the daily minimum tenpereture
  • the daily mean tenperature of 1884-present at the Hong Kong Observatory. The starting dates of the theree data are not same, and so we make use of data of 140 years of 1885-2024. The Hong Kong Observatory dataset has no daily temperature records from 1 Jan 1940 to 31 Dec 1946, due to the disruption of WWII and the Japanese occupation. This 7-year gap means that analyses requiring continuous daily data over multiple years will be affected if these years are included. To ensure consistency, we restrict the analysis period to 1 Jan 1947 – 31 Dec 2024.

Run the following cell to donwload the public data and show the table.¶

In [13]:
Downloads.download("https://data.weather.gov.hk/weatherAPI/opendata/opendata.php?dataType=CLMMAXT&rformat=csv&station=HKO","HK_Maximun_Temperature.csv");
Downloads.download("https://data.weather.gov.hk/weatherAPI/opendata/opendata.php?dataType=CLMMINT&rformat=csv&station=HKO","HK_Minimum_Temperature.csv");
Downloads.download("https://data.weather.gov.hk/weatherAPI/opendata/opendata.php?dataType=CLMTEMP&rformat=csv&station=HKO","HK_Mean_Temperature.csv");

Amax = CSV.read("HK_Maximun_Temperature.csv", DataFrame, 
			 header=["Year", "Month", "Day", "Maximum", "C"],
			 skipto=370,
			 footerskip=3,
			 types=Dict("Maximum" => Float64),
			 drop=["C"],
			 missingstring="***",
			 normalizenames=true); 
delete!(Amax, 48579:nrow(Amax));
Amean = CSV.read("HK_Mean_Temperature.csv", DataFrame, 
			 header=["Year", "Month", "Day", "Average", "C"],
			 skipto=310,
			 footerskip=3,
			 types=Dict("Average" => Float64),
			 drop=["C"],
			 missingstring="***",
			 normalizenames=true); 
delete!(Amean, 48579:nrow(Amean));
Amin = CSV.read("HK_Minimum_Temperature.csv", DataFrame, 
			 header=["Year", "Month", "Day", "Minimum", "C"],
			 skipto=370,
			 footerskip=3,
			 types=Dict("Minimum" => Float64),
			 drop=["C"],
			 missingstring="***",
			 normalizenames=true); 
delete!(Amin, 48579:nrow(Amin));

A = outerjoin(Amax, Amean, Amin,on = [:Year, :Month, :Day])
Out[13]:
48578×6 DataFrame
48553 rows omitted
RowYearMonthDayMaximumAverageMinimum
Int64Int64Int64Float64?Float64?Float64?
118851117.815.313.9
218851220.516.714.2
318851316.615.414.8
418851418.616.114.6
518851520.817.515.1
618851621.118.316.8
718851719.917.915.4
818851817.315.614.2
918851916.915.714.7
10188511018.116.115.1
11188511119.917.215.3
12188511217.916.014.9
13188511317.114.411.9
⋮⋮⋮⋮⋮⋮⋮
485672024122017.714.911.9
485682024122120.216.913.9
485692024122218.015.813.5
485702024122317.516.515.1
485712024122419.117.415.6
485722024122520.618.516.6
485732024122622.920.118.0
485742024122720.919.218.1
485752024122818.816.915.1
485762024122917.415.413.3
485772024123020.417.714.3
485782024123122.619.817.6

Run the following cell to draw the line graphs of Annual maximum, average, and minimum temperatures.¶

Over the long term, temperatures are rising slowly.

In [14]:
A_1947_2024 = filter(:Year => y -> 1947 <= y <= 2024, A);
	Y = combine(groupby(A_1947_2024, :Year),
    :Maximum => maximum => :Ymax,
    :Average => mean     => :Yavg,
    :Minimum => minimum  => :Ymin
		)

	plot(Matrix(Y)[:,2:4],
         size=(500,350),
		 title="Annual max, ave, and min temperatures at HK Observatory",    
		 titlefontsize=9,
		 titlefontcolor="blue",
		 xlim=(-6,84),
		 xticks=([4,14,24,34,44,54,64,74],
        [1950 1960 1970 1980 1990 2000 2010 2020]),
		 ylim=(0,40),
		 grid=false,
		 xlabel="year",
		 ylabel="C",
		 legend=false)
Out[14]:
No description has been provided for this image

Run the following cell to create a movie showing the monthy maximum, average, and minmum temperatures in each year between 1947 and 2024.¶

In [15]:
animhko = @animate for k=1947:2024

M = combine(groupby(A_1947_2024, [:Year, :Month]),
:Maximum => maximum => :Mmax,
:Average => mean    => :Mavg,
:Minimum => minimum => :Mmin)
	
Myy = sort(M[M.Year .== k, :], :Month) 
	
plot(Matrix(Myy[:,3:5]),
    size=(500,350),
    title="Monthly max, ave, and min temperatures at HK Observatory in $k",
	titlefontsize=8,
	titlefontcolor="blue",
	xlim=(0,13),
	xticks=([1,2,3,4,5,6,7,8,9,10,11,12]),
	ylim=(0,40),
	grid=false,
	xlabel="month",
	ylabel="C",
	legend=false)

end

gif(animhko, "HKO_animation.gif", fps = 2)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/HKO_animation.gif
Out[15]:
No description has been provided for this image

7. Central limit theorem for the binomial distribution $B(n,p)$¶

Let $n$ be a positive integer, and let $p$ be a constant satisfying $0 < p < 1$. Denote by $B(n,p)$ the the binomial distribution of $n$ iid experiments $X_1,\dotsc,X_n$ with probability $p$, that is, a discrete probability distribution on $\{0,1,\dotsc,n\}$ such that

$$\operatorname{P}(S_n=k)=\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}, \quad k=0,1,\dotsc,n,$$

where $S_n:=X_1+\dotsb+X_n$.

It is well-known that for any $a,b\in\mathbb{R}$ with $a<b$,

$$\operatorname{P} \left( a \leqq \frac{S_n-pn}{\sqrt{np(1-p)}} \leqq b \right) \rightarrow \frac{1}{\sqrt{2\pi}} \int_a^be^{-x^2/2}dx \quad (n \rightarrow \infty),$$

which is a typical example of the central limit theorem.

Choose $p\in(0,1)$ and run the following cell to create a movie of the central limit theorem¶

In [29]:
# Choose 0<p<1
p=0.3;

x1 = range(-3, 3, length = 1201);
f1 = exp.(-x1.^2/2)/sqrt(2*pi);

animclt = @animate for m=1:150
n1=10*m
z2=zeros(n1+2*max(100,n1));
for l=0:n1
	z2[l+max(100,n1)+1]=Float64(sqrt(n1*p*(1-p))*factorial(BigInt(n1))*p^l*(1-p)^(n1-l)/factorial(BigInt(l))/factorial(BigInt(n1-l)));
end
	
z1 = zeros(length(x1));
c1 = Int64((length(x1)-1)/2);
for k=1:c1
	z1[c1+k]=z2[max(100,n1)+2+Int64(floor(n1*p+(k-1)/200*sqrt(n1*p*(1-p))))];
	z1[c1+1-k]=z2[max(100,n1)+1+Int64(floor(n1*p-(k-1)/200*sqrt(n1*p*(1-p))))];
end

plot(f1,
	title=(L"B(n,%$p) \rightarrow N(0,1), n=%$n1"),
    size=(500,350),
	grid=false,
	ylim=(0,0.45), 
	xticks=([201 401 601 801 1001;],[-2,-1,0,1,2]),
	yticks=([0 1/sqrt(2*pi);],[0,"1/√2π"]),
	xlabel=(L"x"),
	ylabel=(L"y"), 
	edgecolor=:false, 
	legend=:false,
	lw=2,
	color=:blue)
plot!(z1, seriestype=:shape, color=:yellow, alpha=0.6)

end

gif(animclt, "CLT_animation.gif", fps = 10)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/CLT_animation.gif
Out[29]:
No description has been provided for this image

8. Taylor series of $\sin{x}$¶

It is well-known that $\sin{x}$ ($x\in\mathbb{R}$) equals a power series of $s$:

$$\sin{x}=\sum_{k=0}^\infty\frac{(-1)^kx^{2k+1}}{(2k+1)!}, \quad x\in\mathbb{R}.$$

The right hand side is said to be the Taylor series of $\sin{x}$ or the Taylor expansion of $\sin{x}$. More precisely, for any $R>0$,

$$\max_{x\in[-R,R]} \left\lvert \sin{x} - \sum_{k=0}^K \frac{(-1)^kx^{2k+1}}{(2k+1)!} \right\rvert \rightarrow 0 \quad (K\rightarrow\infty).$$

We now observe this convergence. Set

$$S_K(x) := \sum_{k=0}^K \frac{(-1)^kx^{2k+1}}{(2k+1)!}.$$

Run the following cell to create a movie.¶

In [17]:
x = range(-3*pi, 3*pi, length = 301);
f3=sin.(x);
KK=14;
s3=zeros(KK+1,301);
for k=0:KK
	for l=1:301
        s3[k+1,l]=(-1)^k*x[l]^(2*k+1)/factorial(big(2*k+1));
    end
end
	
S3=zeros(KK+1,301);
S3[1,:]=s3[1,:];
for l=1:301
    for k=2:KK+1
		S3[k,l]=S3[k-1,l]+s3[k,l];
	end
end

animtaylor = @animate for K=0:KK
h3=S3[K+1,:];
plot([f3,h3],grid=false,linewidth=2,ylim=(-1.2,1.2),
    size=(500,350),
	title=(L"\sin(x), S_K(x), \quad K=%$K"),
	xticks = ([0 50 100 151 201 251 301;], ["-3π","-2π","-π","0","π","2π","3π"]),
	xlabel=L"x",
	label=[L"\sin(x)" L"S_K(x)"],
	legend=:topleft,legendfont=font(10))

end

gif(animtaylor, "Taylor_animation.gif", fps = 1)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/Taylor_animation.gif
Out[17]:
No description has been provided for this image

9. Fourier series of a step function and a sawtooth function¶

Let $g(x)$ and $h(x)$ be a $1$-periodic step function and a $1$ periodic sawtooth function, that is,

$$\begin{aligned} g(x) & := \begin{cases} 1, &\ x\in[0,1/2], \\ 0, &\ x\in(1/2,1), \end{cases} \qquad g(x):=g(x-n),\ x\in[n,n+1),\ n\in\mathbb{Z}, \\ h(x) & :=x, \qquad h(x):=h(x-n),\ x\in[n,n+1),\ n\in\mathbb{Z}. \end{aligned}$$

Their Fourier series are

$$\begin{aligned} g(x) & \sim \frac{1}{2} + \sum_{k=1}^\infty \frac{2}{(2k-1)\pi} \sin\bigl(2(2k-1)\pi x\bigr), \\ h(x) & \sim \frac{1}{2} - \sum_{n=1}^\infty \frac{1}{n\pi} \sin(2\pi nx), \end{aligned}$$

Set

$$\begin{aligned} S_{2K-1}[g](x) & = \frac{1}{2} + \sum_{k=1}^K \frac{2}{(2k-1)\pi} \sin\bigl(2(2k-1)\pi x\bigr) \quad (K=1,2,3.\dotsc), \\ S_{-1}[g](x) & :=\frac{1}{2}, \\ S_N[h](x) & = \frac{1}{2} - \sum_{n=1}^N \frac{1}{n\pi} \sin(2\pi nx) \qquad (N=1,2,3,\dotsc), \\ S_0[h](x) & := \frac{1}{2}. \end{aligned}$$

On one hand, at the discontinuous points

\begin{aligned} S_{2K-1}[g](k) = \frac{1}{2} & \rightarrow \frac{1}{2} = \frac{g(k-0)+g(k+0)}{2} \quad (K\rightarrow\infty), \\ S_{2K-1}[g](k+1/2) = \frac{1}{2} & \rightarrow \frac{1}{2} = \frac{g(k+1/2-0)+g(k+1/2+0)}{2} \quad (K\rightarrow\infty), \\ S_N[h](k) = \frac{1}{2} & \rightarrow \frac{1}{2} = \frac{h(k-0)+h(k+0)}{2} \quad (N\rightarrow\infty) \end{aligned}

for all $k\in\mathbb{Z}$. On the other hand, we can prove that for any small $\delta>0$

\begin{aligned} \max_{x\in[\delta,1/2-\delta]\cup[1/2+\delta,1-\delta]} \lvert{S_{2K-1}[g](x)-g(x)}\rvert & \rightarrow 0 \quad (K\rightarrow\infty), \\ \max_{x\in[\delta,1-\delta]} \lvert{S_N[h](x)-h(x)}\rvert & \rightarrow 0 \quad (N\rightarrow\infty). \end{aligned}

The partial sums oscillate violately as $N$ increases near the discontinuous points. This is called the Gibbs phenomenon.

Run the following cell to create $g(x)$, $h(x)$, $S_{2K-1}[g](x)$ and $S_N[h](x)$.¶

In [18]:
# Set Maximun number of terms M2
M2=50;

# g
g = zeros(481)
for l=41:241
    g[l] = 1;
end
for l=441:481
	g[l]=1;
end

# S_{2K-1}[g]
x2 = range(-0.1, 1.1, length = 481);
S2=ones(M2+1,481)/2;
l=10
for k=2:M2+1
	for l=1:481
	    S2[k,l]=S2[k-1,l]+2*sin(2*(2*k-3)*pi*x2[l])/(2*k-3)/pi;
	end
end

#h
h = zeros(481);
for l=1:481
	h[l]=x2[l]-Float64(floor(x2[l]));
end

# S_N[h]
S3=ones(M2+1,481)/2;
for k=2:M2+1
	for l=1:481
	    S3[k,l]=S3[k-1,l]-sin(2*(k-1)*pi*x2[l])/(k-1)/pi;
	end
end

Run the following cell to create a movie.¶

In [19]:
animgibbs = @animate for K2=1:M2+1

N2=K2-1

G1=plot(x2,S2[K2,:],
        grid=false,
		linewidth=2,
		xlim=(-0.15,1.15),
		ylim=(-1/8,9/8),
		title="Step Function and its Fourier Series, K=$N2",
        titlefontsize=9,
		xticks = ([0 1/2 1;], [0,1/2,1]),
		yticks = ([0 1/2 1;], [0,1/2,1]),
		xlabel=L"x",
		label=L"S_{2K-1}[g](x)",
		legend=:topright,
		legendfont=font(10))
	plot!(x2,g,
		linewidth=2,
	    label=L"g(x)")

G2=plot(x2,S3[K2,:],
		grid=false,
		linewidth=2,
		xlim=(-0.15,1.15),
		ylim=(-1/8,9/8),
		title="Sawtooth Function and its Fourier Series, N=$N2",
        titlefontsize=9,
		xticks = ([0 1/2 1;], [0,1/2,1]),
		yticks = ([0 1/2 1;], [0,1/2,1]),
		xlabel=L"x",
		label=L"S_N[h](x)",
	    legend=:top,
		legendfont=font(10))
	plot!(x2,h,
		linewidth=2,
	    label=L"h(x)")

plot(G1,G2, layout=(1,2), size=(700,300))

end

gif(animgibbs, "Gibbs_animation.gif", fps = 5)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/Gibbs_animation.gif
Out[19]:
No description has been provided for this image

10. Weyl's equidistribution theorem¶

Let $\lfloor{x}\rfloor$ be the floor function for $x\in\mathbb{R}$, that is, $$\lfloor{x}\rfloor:=\max\{m \in \mathbb{Z} : m \leqq x\},$$ and set

$$\langle{x}\rangle:=x-\lfloor{x}\rfloor.$$

Then $0\leqq\langle{x}\rangle<1$ and $x\equiv\langle{x}\rangle$ mod $1$.

We say that a sequence $\{a_n\}_{n=1}^\infty \subset [0,1)$ is equidistributed in $[0,1)$ if

$$\frac{\sharp\{n=1,\dotsc,N : a<a_n<b\}}{N} \rightarrow b-a \quad (N \rightarrow \infty)$$

Weyl's equidistribution theorem: For any $\gamma\in\mathbb{R}\setminus\mathbb{Q}$, $\{\langle{n\gamma}\rangle\}_{n=1}^\infty$ is equidistributed in $[0,1)$.

We observe Weyl's equidistribution theorem using a bijection

$$[0,1)\ni\langle{x}\rangle \mapsto e^{2\pi{i}\langle{x}\rangle} = e^{2\pi{i}x} = \cos(2\pi{x})+i\sin(2\pi{x}) \in \mathbb{S}^1:=\{z\in\mathbb{C} : \lvert{z}\rvert=1\}.$$

Choose $\gamma^2$ (positive integers, fractions and decimals), and run the following cell to create a movie.¶

In [20]:
# Choose g2=gamma^2
g2=10;

# Upper bound of N
Nmax=250

# Creating sequence
t10 = range(0, 1, length = 101);
c101=cos.(2*pi*t10);
c102=sin.(2*pi*t10);

γ=sqrt(g2);
x10=zeros(Nmax);
y10=zeros(Nmax);

for n=1:Nmax
	x10[n]=cos(2*pi*γ*n);
	y10[n]=sin(2*pi*γ*n);
end

animweyl = @animate for N10=1:Nmax

scatter(x10[1:N10],y10[1:N10],
		title="Weyl's Equidistribution Theorem",
        grid=false,
        xlim=(-1.1,1.1),
        ylim=(-1.1,1.1),
        aspect_ratio=1.0,
        xticks=false,
        yticks=false,
        label=false,
        xaxis=false,
        yaxis=false,
        annotations=[(0,0.3,L"γ=\sqrt{%$g2}",20),(0,0,L"e^{2πiγn}, n=1,...,%$N10",20)],   
	    annotationguide=:auto)
plot!(c101,c102,label=false,linecolor=:magenta)

end

gif(animweyl, "Weyl_animation.gif", fps = 3)	

[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/Weyl_animation.gif
Out[20]:
No description has been provided for this image
In [ ]: